function [vadist,output] = emqnalgo(vadist,Y,H,Q,varargin)

emopt = struct('Display','off','MaxIter',1e6,'FunctionTolerance',0 ...
    ,'StepTolerance',1e-6,'RelativeGradient',1e-6);
while ~isempty(varargin)
    switch varargin{1}
        case fieldnames(emopt)'
            emopt.(varargin{1}) = varargin{2};
        otherwise
            error(['Unexpected option: ' varargin{1}])
    end
    varargin(1:2) = [];
end
outmsg = '...';
rp = nan;
rg = nan;
tic
ll = -inf;
for iter0 = 1:emopt.MaxIter        
    oll = ll;
    [vadist,ll] = emstep(vadist,Y,H,Q);
    dll = ll - oll;
    if strcmp(emopt.Display,'iter')
        fprintf('%d \t %8.12f \t %8.12f  \t %0.0f \n',iter0,ll,dll,toc)
    end
    if dll<0
        error('Algorithm produced a decrease in log likelihood')
    elseif dll<0.5
        break
    end
end

[numvars,numcomp] = size(vadist.mean);
p2p = @(x) vecparm_unc(x,numvars,numcomp);

parms = p2p(vadist);
[newvadist,ll,g] = emstep(p2p(parms),Y,H,Q);
gt = p2p(newvadist) - parms;
S = zeros(size(parms,1));

for iter = (iter0+1):emopt.MaxIter

    op = vecparm(p2p(parms),numvars,numcomp);
    
    oll = ll;

    d = gt - S*g; %typo in Jamshidian and Jennrich (1997) wrong sign;

    a = linesearch(parms,d,p2p,Y,H,Q);
    if isnan(a)
        dp = gt;
        S = zeros(size(parms,1));
    else
        dp = a*d;
    end     

    [newvadist,ll,newg] = emstep(p2p(parms+dp),Y,H,Q);
    dg = newg - g;
    dgt = (p2p(newvadist) - (parms+dp)) - gt;
    dps = -dgt + S*dg;
    dS = (1 + (dg'*dps)/(dg'*dp))*((dp*dp')./(dg'*dp)) ...
        - (dps*dp' +(dps*dp')')./(dg'*dp);
    
    parms = parms + dp;
    g = g + dg;
    gt = gt + dgt;
    S = S + dS;

    dll = ll - oll;

    rg = norm(abs(g).*(max(abs(parms),1)./max(abs(ll),1)),'inf');
    rp = norm(abs(vecparm(p2p(parms),numvars,numcomp) - op)./max(abs(op),1),'inf');
    rll = abs(dll)/max(abs(ll),1);

    if strcmp(emopt.Display,'iter')
        fprintf('%d',iter)
        fprintf('\t %8.12f \t %8.12f',ll,dll)
        fprintf('\t %8.10f',rg)
        fprintf('\t %8.10f',rp)
        fprintf('\t %8.2f',a)
        fprintf('\t %0.0f',toc)
        fprintf('\n')
    end
        
    if (rg<emopt.RelativeGradient) || (rp<emopt.StepTolerance) || (rll<emopt.FunctionTolerance)
        outmsg = 'Stopping criteria met';
        break
    elseif dll<0
        error('Algorithm produced a decrease in log likelihood')
    elseif dll<(eps(ll)^(2/3))
        outmsg = 'Change in log-likelihood less than eps^(2/3)';
        warning('EM terminated. Change in log-likelihood is less than eps^(2/3)')
        break
    end
                      
end   

vadist = p2p(parms);
numparms = size(parms,1);
output = struct('loglike',ll,'iterations',iter,'numparms',numparms ...
    ,'stepsize',rp,'firstorderopt',rg,'message',outmsg);

end
